Temperature dependence of a vortex in a superfluid Fermi gas 



S. Simonucci, P. Picri, and G. C. Strinati 
Division of Physics, School of Science and Technology 
Universita di Camerino, 62032 Camerino (MC), Italy 
and 

IN FN, Sezione di Perugia, 06123 Perugia (PG), Italy 
(Dated: March 22, 2013) 

The temperature dependence of an isolated quantum vortex, embedded in an otherwise homo- 
geneous fermionic superfluid of inflnite extent, is determined via the Bogoliubov-de Gennes (BdG) 
equations across the BCS-BEC crossover. Emphasis is given to the BCS side of this crossover, 
where it is physically relevant to extend this study up to the critical temperature for the loss of the 
superfluid phase, such that the size of the vortex increases without bound. To this end, two novel 
techniques are introduced. The flrst one solves the BdG equations with "free boundary conditions" , 
which allows one to determine with high accuracy how the vortex profile matches its asymptotic 
value at a large distance from the center, thus avoiding a common practice of constraining the vortex 
in a cylinder with infinite walls. The second one improves on the regularization procedure of the 
self-consistent gap equation when the inter-particle interaction is of the contact type, and permits 
to considerably reduce the time needed for its numerical integration, by drawing elements from the 
derivation of the Gross-Pitaevskii equation for composite bosons starting from the BdG equations. 
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I. INTRODUCTION 

Vortices are at the essence of superfluidity and of its 
deep connection with Bose-Einstein condensation (BEG) 
They have thus received considerable interest with 
the raise of ultra-cold dilute trapped Bose gases [2} , where 
they can be generated by setting the trap into rotation Q 
and have been the subject of experimental investigation 
In this context, isolated vortices or even vortex arrays 
have mainly been studied theoretically in terms of the 
Gross-Pitaevskii (GP) equation for the wave function of 
the condensate [1, 0| , which was specifically introduced 
to describe an isolated vortex in an otherwise uniform 
dilute Bose-Einstein condensate. 

Subsequent interest in ultra-cold dilute trapped Fermi 
gases and in the associated BCS-BEC crossover 0, Q 
(whereby a continuos evolution is achieved from a BCS- 
like situation with highly overlapping Cooper pairs, to a 
BEC-like situation where composite bosons form out of 
fermion pairs and condense at sufficiently low tempera- 
ture) has raised the issue of the description of vortices 
in Fermi systems, for which the Pauli principle requires 
one to consider in general a whole set of one-particle 
wave functions instead of a single condensate wave func- 
tion. In this context, isolated vortices (or even vortex 
arrays) have been studied theoretically in terms of the 
Bogoliubov-de Gennes (BdG) equations which were 
introduced as an extension of the BCS approach [l3| to 
describe a non-uniform Fermi superfluid. Experimen- 
tally, arrays of vortices have been detected throughout 
the BGS-BEG crossover once trapped Fermi atoms were 
set into rotation 



involved and time consuming than the solution of the 
GP equation for the bosonic condensate wave function. 
For this reason, consideration has essentially been lim- 
ited to the study of an isolated vortex (with the excep- 
tion of arrays of vortices in the weak-coupling (BCS) limit 
12, 13). In particular, an isolated vortex was considered 
by solving the BdG equations in Refs. 14 1 and [l3| at zero 



From the computational side, solution of the BdG 
equations for the fermionic wave functions is much more 



tem per ature throughout the BCS-BEC crossover, and in 
Ref. jl6l| at finite temperature but in the weak-coupling 
(BCS) limit only. In these works, the superfluid was en- 
closed in a cylinder of radius R. 

Aim of the present paper is to extend the calculation 
of the fermionic BdG equations for a single vortex over 
the whole temperature range from zero up to the critical 
temperature Tc for the loss of the superfluid phase, while 
spanning at the same time the entire BCS-BEC crossover. 
In practice, the crossover between the BCS and BEG 
regimes is essentially exhausted within a range w 1 about 
the unitary limit at {kpap)~^ = where the scattering 
length Of of the two-fermion problem diverges [kp being 
the Fermi wave vector related to the bulk density rip via 
no = kl/iTi'^). 

This will require us to avoid constraining the superfluid 
within a cylinder of radius R with rigid walls, but to let 
it be free of expanding its size without bound when ap- 
proaching Tc from below. To this end, appropriate "free 
boundary conditions" will have to be implemented for the 
BdG equations, in order to recover their correct asymp- 
totic solution far away from the center of the vortex when 
its size would exceed any reasonable value one could take 
for R. The advantage of avoiding the use of a finite value 
R can be perceived, in practice, even somewhat away 
from Tc, as it can be seen from the weak-coupling case 
reported in Fig[l]for the sake of example. 
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FIG. 1. Profile of tfie order parameter A(p) (normalized to 
its asymptotic value Ao at the given temperature) for an iso- 
lated vortex versus the distance p from the center. The weak- 
coupling case with (kpCLF)^^ = —1 is considered for three dif- 
ferent temperatures: (a) T = 0; (b) T = O.eTc; (c) T = 0.9Tc. 
In the three cases, the calculation using "free boundary condi- 
tions" (full line) is compared against that using a cylindrical 
box (dashed line), with the value R = 25fc^^ for the radius as 
typically taken in previous calculations [14| . 



In this way, we will be able to obtain the healing length 
for an isolated vortex as a function of the temperature T 
from T = up to (quite close to) Tc and of the coupling 
parameter {kpaF)~^ spanning the BCS-BEC crossover. 

This information about the way the superfluid healing 
length can be fine-tuned in a Fermi gas, by varying not 
only the temperature but also the inter-particle coupling 
(or both), may also be relevant for the emerging field 
of superfluid interferometers [l3|, in the case it could be 
possible to realize them in practice by coupling systems 
of ultra-cold dilute trapped Fermi atoms. 

In the course of the present calculation, we shall 
also improve on the regularization procedure of the self- 
consistent gap equation which was used in the literature 
for similar problems [l8l - [20j and is required when, like 
in the present context, the inter-particle interaction is of 
the contact type. This will permits us to reduce consid- 
erably the computational time needed for the numerical 
integration of the BdG equations, while leaving unaltered 
the numerical accuracy. To this end, elements will be 



drawn from the derivation of the GP equation for com- 
posite bosons that was obtained in Ref. j2l| on the EEC 
side of the crossover starting from the BdG equations. 

A second, yet not less important, purpose of the 
present paper is to obtain as accurate as possible nu- 
merical solutions of the BdG equations for a non-trivial 
physical problem (like that of an isolated vortex embed- 
ded in an infinite superfluid) under a wide variety of cir- 
cumstances. These numerical solutions could, in fact, 
be used in the future as a "benchmark" for the results 
obtained alternatively by solving approximate local (dif- 
ferential) equations for the gap parameter, which could 
take the place of the fermionic BdG equations at least 
in some approximate sense. In turn, these local equa- 
tions should be better suited to deal with more complex 
problems like the arrays of vortices and the moment of 
inertia of the superfluid, which can be explored exper- 
imentally with ultra-cold trapped Fermi atoms [ll|j [22| 
but remain too difficult to be approached theoretically 
by solving directly the fermionic BdG equations. 

As an example, the GP equation for composite bosons 
at low temperature that was derived in Ref.[2l[ has al- 
ready been tested in the context of the Josephson effect 
psj . to produce results quite similar to those obtained 
by solving the BdG equations on the BEG side of the 
crossover, albeit in a much more efficient way. Sim- 
ilarly, the Ginzburg-Landau (GL) equation for Cooper 
pairs, that was derived by Gorkov [2J| on the BCS side 
of the crossover and close to Tc also starting from the 
fermionic BdG equations, can be most readily applied to 
non-uniform superconductors under a variety of circum- 
stances [2^ since its solution is considerably simpler than 
that of the original BdG equations. 

Along these lines, attempts have already been made in 
the past to derive from the BdG equations extensions of 
the GL equation, which would still apply to the weak- 
coupling (BCS) regime but somewhat dee per in the su- 



perfiuid phase away from the vicinity to Tc More 
recently, a systematic expansion of the BdG equations in 
terms of the small parameter (Tc — T) /Tc was considered 
in the weak-coupling regime, but was explicitly tested for 
the spatially uniform case only [2^. A satisfactory test 
of the above (as well of other) proposal for differential 
equations, that aim at extending the validity GL equa- 
tion deep in the superfluid region, is thus apparently still 
pending and the accurate solution of the BdG equations 
we obtain in the present paper may provide the awaited 
ground for this comparison. 

In this context, an additional important information 
that can be obtained by the present approach comes from 
the analysis of how alternative energy ranges in the so- 
lution of the BdG equations (namely, bound states, and 
near and far continuum) contribute to the different spa- 
tial regions in which the profiles of physical quantities as- 
sociated with a vortex (like the gap parameter itself and 
the number and current densities) can be partitioned. 
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This kind of information is, in fact, expected to be rel- 
evant in future work in order to assess the vahdity of 
approximate local (differential) equations for the gap pa- 
rameter. 

The paper is organized as follows. Section II considers 
the solution of the BdG equations for an isolated vor- 
tex embedded in an infinite superfluid, for which "free 
boundary conditions" are introduced and the associated 
normalization of the wave functions in the continuum is 
obtained. The spatial profiles of the vortex obtained in 
this way under a variety of circumstances are reported in 
Section III. Section IV discusses the procedure through 
which the healing length of the vortex, as a function of 
coupling and temperature, can be extracted from the 
above profiles. Section V provides an analysis of the 
contribution of the different energy ranges in the BdG 
equations to different portions in the profiles of physical 
quantities. Section VI gives our conclusions. The way 
the boundary conditions are implemented is discussed in 
detail in Appendix A, the improved regularization pro- 
cedure for the gap equation is derived in Appendix B, 
and the related expressions for the number and current 
densities are reported in Appendix C. 



II. SOLUTION OF THE BOGOLIUBOV-DE 
GENNES EQUATIONS WITH FREE BOUNDARY 
CONDITIONS 

In this Section, we discuss in detail the solution of the 
fcrmionic BdG equations in cylindrical coordinates for an 
isolated vortex embedded in an otherwise infinite super- 
fluid. To be able to deal with situations when the size of 
the vortex grows without bound upon approaching Tc (in 
practice, when it exceeds a few dozens times kp^), an ex- 
plicit numerical integration of the BdG equations will be 
performed from the center of the vortex outwards only in 
a limited radial range, at the boundary of which connec- 
tion with asymptotic solutions will be sought in terms of 
known functions of mathematical physics. Knowledge of 
these asymptotic solutions will also enable us to to deter- 
mine the normalization of the eigen-solutions of the con- 
tinuum part of the spectrum of the BdG equations. This 
step is of particular importance, since it turns out that 
the continuum part of the spectrum exhausts in practice 
most part of the contribution to the relevant physical 
quantities. 

A. BdG equations for an isolated vortex 
embedded in an infinite medium 



The fermionic BdG equations read: 



H{r) A(r) 
A(r)* -H{r) 



(1) 



where 'H(r) = -V^ /2m — fjL (m being the fermion mass, /i 
the chemical potential, and Ti = 1 troughout). The local 
gap parameter A(r) is determined via the self-consistent 
condition: 

A(r) = -«o W«-(r)* [1 - 2Jp{e,)] (2) 



where /F(e) = (e'^/^'^^"^^ + is the Fermi function at 
temperature T {kB being Boltzmann constant) and vq is 
the (bare) coupling constant of the contact interaction. 
Only positive values of the eigenvalues Ci, can be explic- 
itly considered 0. 

We are specifically interested in a spatially dependent 
gap parameter A(r) with cylindrical symmetry 



A(r) = A(p,vP,2) = A(p)e^' 



(3) 



that corresponds to an isolated vortex directed along the 
z axis with circulation quantum n in integer). [Wc shall 
take n = 1 eventually.] The associated wave functions of 
Eqs.ll]) have the form: 



i^.,£,fejr) = «,(p)e^(^-")^e^'=^- 



(4) 



(£ integer) where A(p). Ui,{p), and Wiy(p) are real func- 
tions. The BdG equations ^ then become: 

Ot Ui,{p) + A(p)i;^(p) = e„u^{p) 

-Ot-n V„{p) + A{p)u„{p) = £yVv{p) (5) 



involving the radial operator 



1 d _ 

2mp dp \ dp 



Imp^ 



(6) 



where fi = p — k'^/2m is the reduced chemical potential. 

Each of the two second-order differential equations ([5]) 
admits a regular solution in p = 0, which behave respec- 
tively as u^{p) ~ p'^' and Vi,{p) ^ pl^^"'. In particular, 
for n ~ 1 (whereby A(p) = r]p for p — > with 77 constant), 
two independent solutions of the coupled equations ([S]) 
can be obtained by taking the indicial conditions: 



,(1) 



«W(p)=/3pl^l 



+3 



where (4 + 3\£\ + 1)13 /m + 77 = 0, and 
^^f(p) = /-^l +••• 



(7) 



(8) 



where - (|£ 3)2]7/(2m) -t- 7/ = 0. 

The differential equations ([S]) are integrated numeri- 
cally from p = up to an outer value i?out, and for 
several values of £ up to a maximum value ^max- Here, 
?out can be related to each other as follows: 
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(i) To begin with, one selects a cutoff energy Ec such that 
only (positive) eigenvalues e^/ up to Ec — arc explic- 
itly considered in the solution of Eqs.® (the remaining 
eigenvalues larger than E^ — fJ- will be dealt with sepa- 
rately by the regularization procedure for the gap equa- 
tion described in Appendix B); 

(ii) One then chooses a value of i?out such that for p > 
i?out the gap A{p) in Eqs.(IS|) has reached its asymptotic 
(bulk) value Aq, say, within 1% (values kpRout — 60-;-200 
prove sufficient for all practical purposes); 

(iii) Finally, one solves numerically Eqs.® for values of 
£ up to ioinx such that i'^^^/ {2mRl^^) ^ E^, that is to 
say, £max kcRout with fc^ ~ \/2mEc (in practice, we 
have taken fmax not smaller than 200). 

It is clear that a reasonable estimate of the value of 
i?out entails knowledge of the profile of A(p), which in 
turn requires the solution of the self-consistent condition 
(HI). Wc defer to Appendix B the solution of Eq. © to- 
gether with a proper treatment of the convergence of the 
sum over v for large values of e^. In this context, a novel 
regularization procedure for the gap equation ^ will be 
introduced, which improves on regularization procedures 
previously considered in the literature 3 201 (thereby 
effectively reducing the numerical value of Ec). 

B. Asymptotic behavior of the wave functions 

For an isolated vortex embedded in an otherwise infi- 
nite superfiuid medium, the eigenvalues ei, of the BdG 
equations ([Sj belong to a continuous spectrum above 
the threshold Aq (apart from the Andreev-Saint- James 
bound states that lie below this threshold) . For the wave 
functions belonging to this continuum, in turn, the nor- 
malization is determined from their "asymptotic" behav- 
ior for large values of p, which may be identified only 
for p ^ Rout- For this reason, the asymptotic behav- 
ior of u^{p) and v^(p) for p — > oo has eventually to be 
searched in terms of known functions of mathematical 
physics, which is however not possible for the radial BdG 
equations ([5]) as they stand. 

To overcome this problem, we have adopted the fol- 
lowing strategy. If i?out is large enough, the centrifugal 
terms £^/(2mp^) and (£ — n)^/(2m/9^) in Eqs.® are im- 
portant only for large values of £, in such a way that we 
may replace i and {£ — n) by their average value: 



_ [t+{£-n)] n 
2 2 



£- 



1 



(9) 



By this replacement, in the centrifugal terms for p > i?out 
we make an error smaller than Ec/tmax- Accordingly, for 
p > Rout in the place of Eqs.(IS]) we consider the following 
"modified" BdG equations with a common value of £': 



Oi' Uy{p) + Aq Vy{p) = El, u„{p) 

-Of Vy{p) + Ao u„{p) = El, v^{p) . 



These coupled equations can be solved analytically in 
terms of known functions of of mathematical physics, by 
considering the auxiliary equation: 



1 



2mp dp \ dp J 2mp 



fk{p) 



2m 



fJ- j fk{p) 
(11) 



where fc^ = fc^^ + fc^. This equation is equivalent to the 
canonical equation of the Bessel functions of index \£'\ 







(12) 



in the dimensionless variable ^ = kj^p [30| . The solutions 
to Eqs. ([T(I|) arc thus sought in the form 



u^{p) = Uk f{k±p) 



.ip)=vkf{k^p) , (13) 



which reduce Eqs.(|TO|) to the standard system of alge- 
braic equations [9| 



p] Uk + AoVk ^ Ek Uk 

2m 



£k Vk , 



p] Vk + Ao Uk 

2m 



yielding 



wf, = ; 1 



= l-vt 



(14) 
(15) 




When dealing with the continuum spectrum, it is con- 
venient to use the energy eigenvalue e as the independent 
variable. This constrains k± in Ea. (|T6)) to the values: 



(17) 



fc^ = ± \ 2mp ± 2m 



E^-Al 



(10) 



for given e and k^ . To comply with the notation originally 
introduced in Ref.jsij to describe tunneling through a 
barrier in a superconductor, wave vectors with the plus 
(minus) sign inside the square root in Eq. (|17p are referred 
to as electron-like (hole-like) wave vectors. 

Depending on the value of e and the sign of p, there can 
be alternatively four complex solutions, four real solu- 
tions, and two real and two complex solutions of Eq. (|17|) . 
Only complex solutions resulting in decaying exponen- 
tials for (0 — > oo can be accepted. A discussion of the 
explicit solutions in the various energy ranges, depending 
also on the sign of fl, is reported in Appendix A, where 
the boundary conditions at p = Rout between the numer- 
ical solutions of Eqs.(IS]) for p < Rout and the analytical 
solutions of Eqs. ([T(I|) for p > i?out are also reported. 
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What is relevant here is that, depending on the allowed 
solutions k± to Ea. ((T7)) . the solutions {u^{p),v^{p)) of 
the BdG equations for p > i?out can be expressed as 
linear combinations of Bessel Jq(C), Neumann ¥^(0, and 
Hankel i?^(C) functions of index a = \£'\ and argument 
C = k±p. These functions, in turn, have the following 
asymptotic behaviors (that holds for p S> -Rout ) [iOl ■ 



cos 



, 1 1 

C — —Tia -TT 

^2 4 

. 1 1 

C 7ra TT 

^2 4 



(18) 



exp 



1 



-7ra - 



The behaviors ([T8|) are what is only needed to calculate 
the normalization of the wave functions in the continuum 
part of the spectrum, to be considered next. 

C. Normalization in the continuum 

The normalization of the (two-component) wave func- 
tions, that are solutions of the BdG equations ([5]) for en- 
ergy eigenvalues lying in the continuum, can be obtained 
by ada ptin g to the present context the method discussed 
in Ref.[32| for the Schrodinger equation. 

Let us consider the BdG equations ([S]) for two different 
energies e and e', both lying in the continuum. Multi- 
plying these equations from the left by the pair (it^/ , v^' ) 
and (itejUe), in the order, subtracting the resulting ex- 
pressions side by side, and integrating over the radial 
coordinate from p = up to p = i?. we obtain: 



dpp u^,{p)u^{p) 



R 



2m (£-£') 



■ip) 



dulip) 



(19) 



dp 
dp 



<{P) 



Veip) 



du^, jp) 
dp 

dv^',{p) 



dp 



where the index A distinguishes degenerate independent 
solutions (cf. Appendix A) and an integration by parts 
has been performed. In the expression the two 

limits i? — > c» and e — e' have been taken in the or- 
der. Note how the division by (e — e') is interpreted as a 
principal part value (7^), consistently with the "standing- 
wave boundary conditions" we are adopting for the ra- 
dial problem. Note further that the (extreme) asymp- 
totic form of the wave functions is what is only needed 
to establish their normalization. 

In particular, for an asymptotic form of the type (with 
real values of fcj^): 



v^{p) 



Uk 



[cMk±p) + d^Y^{k±p)] (20) 



where c\ and d\ are real coefficients, in the appropriate 
limits the expression ((T9)) reduces to: 



dpp 



u^,ip)u^ip) + v^,{p)v^ip) 



[c\c\' + dxdy] — S{k±_ - k'^) . 



(21) 



To obtain this result we have made use of the identity: 



1 



£ (e-e') 



1 



k^ {k^-k'^) 



(22) 



that holds in the limit e — > e'. A simple generalization of 
the expression (j2ip can be obtained when more than one 
wave vector appear on the right-hand side of Eq. (PU)) . 

III. SPATIAL PROFILES OF A VORTEX FROM 
ZERO TO THE CRITICAL TEMPERATURE 

The solution of the BdG equations for an isolated vor- 
tex embedded in an infinite superfluid, discussed in Sec- 
tion II, enables us to obtain the spatial profile A(p) of the 
gap parameter (via the regularized gap equation ([50)) of 
Appendix B), as well as of the number n[p) and current 
i (p) densities (whose asymptotic contributions are given 
by Eqs. ([54)) and ([55]) of Appendix G, respectively). 

In the following, the chemical potential p entering the 
BdG equations is eliminated in favor of the asymptotic 
(bulk) value uq of the density via the standard BCS den- 
sity equation for a homogeneous system in the absence 
of the vortex, namely. 



no 



dk 

(2^ 



(23) 



where 



2m 



p and = y^^ + A^, since correc- 



tions to p due to the presence of an isolated vortex are 
negligible in the thermodynamic limit. This procedure, 
in turn, fixes the value of /cf- 

Figures show our numerical results for the quanti- 
ties A(p), n(p), and j{p), in the order, for the four cou- 
plings (kpaFy^ = (-2.0,-1.0,0.0,-hl.O) and the three 
temperatures T = (0.0, 0.5, 0.9)Tc. These plots were gen- 
erated using a common cutoff energy Ec = SEp which, 
thanks to our novel regularization procedure (cf. Appen- 
dices B and G), proves sufficient to achieve maximum ac- 
curacy of the calculations to the extent that using larger 
values of Ec provides essentially the same results. A num- 
ber of similar plots (not shown here) have also been sys- 
tematically generated over a finer mesh of temperatures 
from T = up to r = 0.95Tc, in order to extract from 
them the temperature dependence of the healing length 
associated with the vortex, as discussed in Section IV. 

Note from Figsl2]|4] that the size of the vortex in- 
creases more rapidly with increasing temperature when 
approaching the BGS limit (kpaF)'^ ^ —1. Note also 



6 




FIG. 2. Gap parameter A(p) (normalized to its asymptotic 
value Ao at the given temperature) of an isolated vortex ver- 
sus the distance p from the center, for the coupling (fc_Ftt_F)~^: 
(a) -2.0; (b) -1.0; (c) 0.0; (d) +1.0. For each coupling, 
three different temperatures are considered: T = (full lines); 
T = O.STc (dashed lines); T = O.OTc (dashed-dotted lines). 



the presence of the characteristic Friedel's oscillations 
in all these quantities when this limit is approached at 
low temperature. These oscillations, however, fade away 
rather quickly as the temperature is increased toward Tc- 

As we have already mentioned, the reason why we 
have invested much effort in determining the continuum 
part of the spectrum of the BdG equations in an infinite 
medium is that this part is expected to exhaust in prac- 
tice most part of the contribution to physical quantities. 

In support to this expectation, we show in FigIS] the 
profiles of A(p), n(p), and j{p) obtained at zero tem- 
perature for the coupling {kpap)~^ — —1, alternatively 
by including or omitting the contribution from the con- 
tinuum part of the spectrum in the calculation of these 
quantities. Drastic changes in these profiles result indeed 
when the contribution from the continuum is omitted 
from the calculation (with similar conclusions drawn for 
different temperatures and couplings). A more complete 
analysis of how different energy ranges in the solutions 
of the BdG equations contribute to the spatial profiles of 
these physical quantities will be presented in Section V. 
Note that an appropriate absolute normalization is used 



FIG. 3. Number density n{p) (normalized to its asymptotic 
value no) for an isolated vortex versus the distance p from 
the center, for the same couplings and temperatures of Fig[2] 
The inset shows the density at the center of the vortex vs 
(^fsf)"^, where our results at T = (0, 0.5, 0.9)Tc from bot- 
tom to top (circles with interpolating dashed lines) are com- 
pared with those at T = from Ref.[l4| (squares) and from 
Ref.[li| (triangles). 



in Figl5]for each quantity, in order to obtain a meaningful 
comparison. 

It is further relevant to compare the profiles of the or- 
der parameter and the number current obtained by the 
present accurate solution of the BdG equations on the 
BCS side of the crossover close to with those ob- 
tained by the less demanding numerical solution of the 
Ginzburg-Landau (GL) differential equation for the order 
parameter Aql, namely [33| : 



1 



7C(3)-Bf 



4£; 



T 

T'c 
GL 



Am 



AglW 



(r)=0 



(24) 



where C(3) — 1.202 is the Riemann zeta function of argu- 
ment 3. In terms of this AglIi"), the GL current is then 
given by the expression [33|: 



jGL(r) 



7C(3)no 



16 i TO {TrksTc 



-[AGL(r)*VAGL(r) 
-AglWVAglW*] .(25) 
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FIG. 4. Number current j{p) (normalized to its maximum 
value jmax at the given temperature) of an isolated vortex 
versus the distance p from the center, for the same couplings 
and temperatures of Figs [2] and [S] The maximum value of the 
current conventionally identifies the vortex radius Rv. 



Since the equation (|24| for Aql and the expression ([25|) 
for Jgl have been derived microscopically from the BdG 
equations in the (extreme) BCS limit and close to the 
critical temperature , one expects the numerical com- 
parison with the full solution of the BdG equations to im- 
prove as these limiting conditions are approached. That 
this is indeed the case is shown in FigsIS] and [71 where 
already for the coupling (kpap)^^ = —2 and the temper- 
ature T = 0.95Tc the comparison between the two (BdG 
and GL) calculations appears quite good. 

The above example can be regarded as a prototype for 
what was meant in the Introduction, about the fact that 
non-trivial numerical solutions of the BdG equations can 
be used in practice to test the validity of local equations 
for the order parameter under specific circumstances. 

IV. EXTRACTING THE TEMPERATURE 
DEPENDENCE OF THE COHERENCE LENGTH 

From the spatial profiles A{p) of the gap parameter 
for an isolated vortex that were obtained in Section III, 
we can now extract the characteristic coherence (heal- 
ing) length as a function of temperature and coupling 
according to the following procedure. 



FIG. 5. (a) Gap parameter A(p), (b) number density n{p), 
and (c) number current j(p) of an isolated vortex versus the 
distance p from the center, at zero temperature for the cou- 
pling (kpaF)^^ = —1. The results of the complete calcu- 
lation (full lines) are contrasted with those of a calculation 
that excludes the contribution from the continuum part of 
the spectrum in Eqs. (|50p . (|52p . and (|53|l (dashed lines). 



Wc note at the outset that, for given temperature, 
A(p) approaches its asymptotic (bulk) value Aq far away 
from the center of the vortex with the power-law behav- 
ior Ao(l — C^/2/o^), where ^ is a characteristic length. In 
particular, in the BCS limit close to T^, this behavior can 
be obtained directly from the GL equation (1241) whereby 
C is identified with the GL coherence length [33| : 



?gl(T) 



7C(3) Ef 



1 



24 m nksTc 



T 



-1/2 



(26) 



Similarly, in the BEG limit close to zero temperature, one 
can resort to the GP equation for composite bosons onto 
which the BdG equations map in that limit jU] , and iden- 
tify Q with the GP healing length = (87rapno)^^/^. 

Quite generally, for any coupling and temperature 
smaller than Tc, we have verified from the numerical solu- 
tion of the BdG equations that A(/9) always approaches 
its asymptotic value Ag like p~'^. In practice, we have 
obtained the value of C, through a fit of the type: 
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T=0.95T(, , (kpap) =-1 



T=0.95T(, , (kpap) =-2 





FIG. 6. (a) Gap parameter A(p) and (b) number current j{p) 
of an isolated vortex versus the distance p from the center, 
close to the critical temperature for the coupling (kpaF)'^ = 
— 1. The results of the calculation of the BdG equations (full 
lines) are compared with those obtained by the Ginzburg- 
Landau (GL) theory (dashed lines) . The maximum values Ao 
for A(p) and jmax for j(p) correspond to the BdG calculation. 



A(p) = CO 1 



2p2 



when Ai?v < p < (50-^150)fc^^ 

(27) 



with A ~ 2 6 depending on coupling and temperature. 
Here, i?v is the vortex radius identified from the profile of 
the current like in Fig]?] For smaller values of p, however, 
we have found that a separate exponential fit of the form 



A{p) = 6o 1 - 5i e 



when kp^ < P < Ai?v 



(28) 



is more appropriate. The need to exclude values of 
kpp smaller than one (at least on the BCS side of the 
crossover) , in order to identify the length ^ as in Eq. (^5]) , 
was pointed out in Ref . [l^ for an isolated vortex and in 
Ref.[23| in the context of the Josephson effect. 

The two independent fits ([27| and p8)) determine the 
two length scales C and ^ which may, in principle, be 
different from each other. We actually expect the ra- 
tio ^/C not to be appreciably different from unity for all 
couplings and temperatures, in such a way that a sin- 
gle length scale can be eventually identified also from the 



FIG. 7. Same as Fig[6]but for the coupling {kpap) ^ 
closer to the BCS limit. 



BdG equations. This would be similar to what occurs 
both in the BCS limit close to Tc and in the BEC limit 
close to zero temperature, where a single length scale 
{^gl{T) and ^gp, in the order) is identified. 

Figure [5] shows the typical quality of the fits (P?]) and 
(P5| in the two adjacent spatial regions, for a specific cou- 
pling and three different temperatures. These fits have 
then been repeated for several couplings about unitarity 
and for a rather dense mesh of temperatures. 

The values of the healing length ^ extracted from 
these fits, for several couplings and from T = up to 
T = 0.99Tc, are reported in Fig|9] as black dots. Note 
again how ^ increases faster with increasing temperature 
when the coupling progresses toward the BCS limit. The 
dashed lines in FiglH] are then obtained by assuming a 
simple expression of the form kF^{T) = A{1 — T/Tc)~^^'^ 
to hold for any coupling over the whole temperature 
range from T = up to (very close to) Tc, in analogy 
to the GL expression kp^GhiT) = ^gl(1 - T/Tc^^/^ 
[for which ^GL = 0A7Ef/Ao{T = 0), cf. Eq.^] that 
holds in principle only in the (extreme) BCS limit quite 
close to Tc- From this kind of fit we obtain the values A = 
(13.41, 3.08, 0.96, 0.73) for the four couplings (fc^aF)"^ = 
(—2.0. — 1.0, 0.0, +1.0), in the order, which can be com- 
pared with the GL values ^gl = (10.12,2.26,0.68) for 
the three couplings (fc^af)"^ = (—2.0. — 1.0, 0.0), values 
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FIG. 8. Gap parameter A(p) (normalized to the Fermi en- 
ergy Ef) of an isolated vortex versus the distance p from the 
center, for three different temperature and {kpciF)'^ = — 1- 
The results of the fittings to extract the lengths ^ and in 
the two different intervals of p (broken lines) are compared 
with those of the the full calculation (full lines) . The value of 
the vortex radius Rv is marked in each case. 



FIG. 9. The values of the healing length ^, as obtained from 
the fits (|28p to the profiles of A(p) (dots) for four different 
couplings across unitarity, are shown versus the temperature 
T (in units of the respective critical temperature Tc). These 
values are then fitted by the mean-field-like expression (,{T) oc 
{Tc — T)~^^'^ over the whole temperature range down to T = 
(dashed lines). 



that are determined only in terms of the corresponding 
values of Ao(r = 0). 

Similar plots can be produced for the other length scale 
C extracted from the fits ([?7)) to the profiles of A{p). 
Figure [TOl shows the ratio between these two length 
scales as a function of temperature for several couplings. 
It is rather remarkable that this ratio remains quite close 
to unity in all cases we have considered, thus justifying 
the statement that a single length scale (say, the healing 
length ^ of Ea.ipS])) can meaningfully be extracted from 
the BdG equations for all couplings and temperatures. 
This conclusion will also be confirmed by a similar anal- 
ysis about the temperature and coupling dependence of 
the vortex radius Rv reported in Appendix C. 

Finally, it is interesting to compare the values of the 
healing length ^ at zero temperature across the BCS- 
BEC crossover, obtained by the present BdG analysis 
of the spatial profile of the gap parameter for an iso- 
lated vortex, with the alternative (and, in principle, un- 
related) results for the so-called "phase" coherence length 
^phase, which wcrc obtained in Ref. [s^l from the analy- 



sis of the spatial variation of the longitudinal component 
of the correlation function of the order parameter in an 
otherwise homogeneous system. This comparison, pre- 
sented in Fig llli shows a remarkable overall agreement 
between the coupling dependence of these two quanti- 
ties, for which the minimum occurs at about unitarity 
in both cases. In addition, the inset of Fig[TT] presents 
similar curves obtained by the present BdG analysis at 
finite temperatures. In this case, the minimum is seen to 
move for increasing T progressively toward the BEG side 
of unitarity, as it is expected from the slower increase of 
the healing length ^ for increasing temperature when the 
coupling progresses toward the BEG limit. 

V. CONTRIBUTION TO THE SPATIAL 
PROFILES OF PHYSICAL QUANTITIES FROM 
DIFFERENT BDG ENERGY RANGES 

We have already pointed out in Section III (see FigE] 
therein) that the continuum part of the spectrum of the 
BdG equations contributes in a substantial way to the 
spatial profiles of the gap parameter as well as of the 
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four different couplings across unitarity versus the temper- 
ature T (in units of the respective critical temperature Tc). 
The horizontal (dashed) lines mark the value unity for ^/(^. 



number density and current. 

In particular, we have obtained the result that the 
bound-state part of the spectrum which lies below the 
continuum threshold does not contribute to the density at 
the center of the vortex, so that in this case the contribu- 
tion of the continuum part of the spectrum is overwhelm- 
ing. On the other hand, from the form of the analytic 
result ((54)) for the asymptotic contribution to the density 
that originates from the continuum levels at high energy, 
one concludes that these levels, too, do not contribute to 
the density at the center of the vortex since |A(r)| van- 
ishes therein. It thus appears interesting to determine 
the way different energy ranges in the solutions of the 
BdG equations contribute to the spatial profiles of the 
above physical quantities. 

To this end, wc introduce an "upper limit" E^i for 
the energy such that only eigenstates of the BdG equa- 
tions with El, < i?ui ~ ^J■ are retained in the calculation 
of the partial gap parameter and of the partial num- 
ber density and current. We then increase E^^i progres- 
sively starting from its value at the continuum threshold, 
which corresponds to E^i — /i = Aq when /i > and to 
Eui — M = Aq + /i^ when /i < 0, reaching large values 
of Eui to include eventually the high-energy part of the 
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FIG. 11. Healing length ^ at zero temperature versus the cou- 
pling (kFap)^^ , obtained from the present BdG calculation 
(full line) and from the approach of Ref. [s^] referred to as PS 
(dashed line). For the sake of comparison, the two curves are 
made to coincide in the extreme BCS limit by multiplying the 
BdG curve by a factor 1.257, in order to account for their dif- 
ferent definitions. In addition, the inset shows ^ vs (Afcif)"^ 
obtained from the BdG calculation at the finite temperatures 
T = 0.2Tc and T = O.STc across the BCS-BEC crossover. 



continuum. 

Accordingly, the partial value A< (r) of the gap param- 
eter, that includes only eigenstates up to iJui — ^J■, can be 
obtained from Eg. ([50]) of Appendix B by discarding the 
contribution of eigenstates with energy above E^i — /i, 
thus writing in the place of Eg. ([50)1 : 



u,ir)v,{rr[l-2fF{e,)] (29) 



where TZ{kc) is defined by Eq. (|40| of Appendix B. Corre- 
spondingly, the partial values n^(r) for the density and 
j<(r) for the current are obtained from the expressions 



(f52|) and (fSB]) reported in Appendix C, where now the 
is limited to energies Si, < E^i — /i. For internal consis- 
tency, however, the eigenstates u^{r) and v,y{r) utilized 
in these partial expressions are calculated from the BdG 
equations with the correct self-consistent value of A(r) 
which includes the contribution from all eigenstates. 

The result of this calculation at zero temperature is re- 
ported in FiglT2] for three characteristic couplings. Par- 
ticularly striking appears here the result for the density 
for the couplings —1.0 and 0.0 [cf. panels (d) and (e) 
of FigllJ] , for which the finite value at the center of the 
vortex is mostly contributed by continuum eigenstates 
quite close to threshold while no contribution is provided 
by the bound states below threshold. In particular, when 



-1.0 the value of n^{p = 0) passes from 20% 



to 99% of its full value when Eui varies from l.^OEp to 
2.0£;f; and when {kpap)''^ = 0.0 from 40% to 95% of 
its full value when E^i varies from lAQEp to 3.0Ep. 
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FIG. 12. Partial radial profiles of the gap parameter A^(p), 
density n^(p), and current j^(p) at zero temperature for 
three couplings ikFap)"^ = ( — 1.0, 0.0, +1.0), obtained us- 
ing different values of the upper limit _Eui for the energy. For 
convenience, the correspondence between the various values 
of Eui and the types of lines used in these plots for the three 
different couplings is reported separately in Table |I] below. 



The above finding, that no contribution to n^{p = 0) 
originates from the bound states, is related to the fact 
that all bound states turn out to correspond to the second 
type of solutions ([8|) with ^ < 0, such that all v^\p) (and 
thus the density) vanish when p ~ 0. At the same time, 
these bound states contribute in a coherent fashion to the 
anti-clock- wise circulation of the current (cf. Eq. ([5^ of 
Appendix C), such that their contribution to the current 
may even exceed the value of the total current which 
includes also the contribution from the continuum [cf. 
panels (g) and (h) of FigHJ]. That the contribution of 
the (Andreev) bound states may sometimes exceed 100% 
of the total current was already poi nted out in the context 
of the Josephson effect in Refs. 23l. l35l|. 
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TABLE I. Correspondence between the values of Eui (in units 
of Ef) and the types of lines used in Fie: ll2l for the three 
different couplings there considered. In all cases, the smallest 
value of Eui corresponds to the continuum threshold. 



VI. CONCLUDING REMARKS 

In this paper, considerable efforts have been devoted 
to obtain an accurate numerical solution of the fermionic 
BdG equations for a non-trivial but still manageable 
problem of an isolated vortex embedded in an otherwise 
infinite superfluid. We have spanned the whole BCS- 
BEC crossover as a function of temperature up to Tc, in 
such a way that the spatial extension and the detailed 
shape of the vortex changes considerably as a function of 
both coupling and temperature. To this end, we have left 
the vortex free of expanding out in the bulk of the super- 
fluid not being constrained by walls, and implemented 
for the purpose the use of free boundary conditions for 
the BdG equations. We have also introduced a new reg- 
ularization procedure for the gap equation that improves 
on previous proposals, so as to reduce the computational 
time while leaving unaltered the numerical accuracy. 

In this way, we have obtained the healing length for 
the vortex structure of the gap parameter as a function 
of both the temperature (from T = essentially up to 
T = Tc) and the coupling parameter {kpap)~^ . This 
quantity shows an interesting behavior across the BCS- 
BEC crossover, which generalizes over the whole temper- 
ature vs coupling phase diagram what is already known 
from: (i) The GL approach in the weak-coupling (BCS) 
limit close to T^; (ii) The GP equation in the strong- 
coupling (BEG) limit at zero temperature; (iii) The ap- 
proach of Ref.[3J] across the BGS-BEG crossover at zero 
temperature. 

In addition, by the present approach we have now at 
our disposal an accurate numerical solution of the BdG 
equations obtained for a non-trivial problem under a va- 
riety of circumstances, against which one might be able 
to compare the results of approximate differential equa- 
tions that originate from local approximation of the BdG 
equations themselves. These local (differential) equa- 
tions could be, for instance, of the GL type in the weak- 
coupling (BCS) limit close to Tc [2J| , or of the GP type in 



the strong-coupling (BEG) limit at zero temperature [21 
In particular, still long awaited appears to be the com- 
parison with the results obtained in the weak-coupling 
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(BCS) limit away from Tc deep in the superfluid phase, 
where generalizations of the GL equation have been at- 
tempted 26 -2^ and deviations between the solutions of 
the BdG equations and these local equations are expected 
at low enough temperature. 

The practical advantage of these differential equations 
stems from the fact that they are considerably simpler 
to solve than the original BdG equations, in such a way 
that, once their validity would have been explicitly tested 
against the results of the BdG equations in a number of 
manageable problems, they could be applied with con- 
fidence to the solution of more complex physical prob- 
lems for which the use of the BdG equations remains 
prohibitive. Work along these lines is in progress [36| . 
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APPENDIX A: ENFORCING THE BOUNDARY 
CONDITIONS AT iJout 

In this Appendix, we describe in detail the solutions of 
the form (fT3|) that hold for p > i?out and are associated 
with the alternative values of k± obtained from Eq. ([T7| , 
depending on the value of the energy e and the sign of 
//. The solutions determined in this way for p > i?out 
will then be used to specify completely the wave func- 
tions obtained numerically for p < i?out, by enforcing 
the appropriate boundary conditions ai p ~ i?out- 

The method we use here is similar to that discussed in 
Ref.pst for a one-dimensional geometry appropriate for 
the study of the Josephson effect throughout the BCS- 
BEC crossover, which extends the original approach of 
Ref.jsH that was limited to the (extreme) BCS limit. 
A related approach was also used in Ref.[37[ for a gap 
parameter with spherical symmetry, and then utilized in 
Ref.jsit to obtain the profile of a single vortex without 
enclosing it in a cylinder, again in the (extreme) BCS 
limit with a large coherence length at zero temperature. 

When the positive energy e is increased from zero past 
the value \/ p"^ + Aq, the four solutions for k± given by 
Eg. p?)) move in the complex fcj^-plane. To follow their 
evolution versus e, it is convenient to label these four 
solutions separately by adopting the convention: 



k^l^ = + \j2mp + 2m i/e^ - 



u(2) 



2mp + 2771 J e'^ - Ag 



FIG. 13. Dispersion relation e vs kj_ given by Eq. (|16P . with 
(a) /i > and (b) jl < 0. These plots identify the six energy 
ranges I- VI, where proper selections of the wave vectors (|30[) 
have alternativefy to be done. 



2m/i — 2to y e 



^0 



iTTip 



2m 



(30) 



For p > i?out, these wave vectors enter the arguments of 
the Bessel Ja{k±p), Neumann Ya{k±p), and Hankel func- 
tions H^{k±p) functions (where a = with asymp- 
totic behavior p^ . according to the following scheme. 

For the six ranges that can be identified depending also 
on the sign of p, (as shown in Figs fTST a) andfT^lb). in the 
order), we obtain by inspection of the expressions (j30p : 



Range I: p. > and < £ < Aq. 



The fc'i^^ {i 



(1) 



= l,---,4) are all complex, but only k]_ 
and k^f^ have a positive imaginary part. We then take 



alternatively A;^^'' and fc^'' in the function H^{k±p). 



;.{4) 



Range 11: p> and Aq < £ < ^/P + Af. 

In this case, /c^' {i = 1, • • • , 4) are all real, and we take 

A:^^' and /c^' in both functions Ja{kj_p) and Ya{k±p). 



Range III: p. > and e > \/ p^ + Aq 



Here, fcl^'' and k'f'' 



are real, and k^f^ and kK^' are purely 



(4) 



imaginary with Imj/c^'} > 0. Wc thus take k''^' in both 



(1) 
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functions Ja{k±p) and Ya{k±p), and in the function 



Range IV: /i < and < e < Aq. 

Same as for range I. We thus take alternatively fc^'' and 
/c^'' in the function H^{k±p). 



Range V: fi <0 and Aq < e < + 

In this case, all the k^f' {i = l,---,4) are purely 

imaginary but only Imlfc^^^} and Imj/c'^^} are positive. 

Again, we take alternatively and fc^'' in the function 

Htik^p). 



Range VI: fl <0 and e > ^//^ + A^. 
Here, fc^^'' and k^f^ are real and fc^-* and fc^-* are purely 
imaginary, but only Im{fc'|^''} is positive. Accordingly, we 
take k^^^ in both functions Ja{k±p) and Ya{k±p), and 
k'f^ in the function H^{k±p). 

With these premises, we pass now to enforce the bound- 
ary conditions at p = i?outi between the numerical solu- 
tions of the BdG equations ([5]) for p < Rout discussed in 
subsection II- A and the analytic solutions of the modified 
BdG equations PU)) for p > Rout introduced in subsec- 
tion II-B. 

Ranges I, IV, and V as specified above can be dealt 
with in the same way, by writing the boundary conditions 
in the form (in the following equations, by £' we shall 
actually mean its absolute value \£'\): 



4^' 
„(i) 



Uki 
Vk, 



(-Rout) 
(-Rout) 



+ b 



(-Rout) 



(2) 



(31) 



Vki 



H+{k'^fRont 



for the functions, and 



dp / p=iiout 

\ dH+{k^l'>p) 



+ b 



dp 
dp 



(32) 



Vkt 



dp 



+ d 



-Rout 



u,, ^| d^+(fcf p) 



dp 



for their first derivatives. Here, {ui^\p),vi^\p)) and 
{ui (p), Ve (p)) are the two independent solutions of the 
BdG equations ([5]) identified by the indicial conditions 
(HI) and (HI, in the order. The conditions (jSl]) and (|32|) 
thus provide an algebraic homogeneous system of four 
equations in the four unknowns (a, b, c, d), which admits 
nontrivial solutions only for special values of e, which cor- 
respond to the Andreev-Saint- James bound states asso- 
ciated with the spatial depression of the gap A(p) about 
p = 0. In this case, the normalization of the single wave 
hmction (^^(p), ^^(p)), as obtained by the linear com- 
bination on the left-hand side of Eq. (PT|) for p < -Rout 
and on the right-hand side of Eq. ipTj) for p > -Rout, is 



determined by: 



dpp [Ue{p)Ue{p) + Ve{p)Ve{p)] = 1 



(33) 



Ranges HI and VI can as well be treated on the same 
footing, by writing the boundary conditions in the form: 



U^e'^Rou 



uY'iRout) 
4''(J?out) 



Uki 
Vk, 



(-Rout) 

cJe,{k^l^R^ut) + dYe,{k^l^Raut) 



+ e 



Uk, 
Uki 



H2,{kfR, 



(34) 



and 



+ e 




dp 



dp ^ P=-Bout 



+ d 



dp 



-Ron 



(35) 



where i = 3 in range HI and i = 4 in range VI. The con- 
ditions and correspond to four algebraic equa- 
tions in the five unknowns (a, fe, c, d, e). The normaliza- 
tion condition (j2ip with A = A' for e in the continuum 
then provides a fifth condition for the coefficients c and 
d, that permits to determine all coefficients uniquely. 

Finally, range II requires a slightly different handling 
because the electron-like and hole-like wave vectors are 
both real. We then apply the boundary conditions to 

the functions {u^e\p)-,v'^\p)) and {u^e\p),vf'\p)) sep- 
arately and write: 



4^' 



(-Rout) 
«r'(J?out) 



(36) 



and 



Uki 

Uk^ 
Vk^ 



dp 
dp 



c^iJi,{k^l^Ro^t) + diiYe,{k'^l'>R^^t] 
ci3 Jr(fcV^'-Rout) + rfi3y«'(fci^'iiout] 



(37) 



P— -Rout 

(iJ^-(fc'^V) 



Cl3 



dp 



dp 



P=-Rou 



P= flout 



dli 



d,i3 



dp 

dYi,{kfp) 



dp 



P = -Rout 



p=-Ro 



Here, the four coefficients (cn, dn, C13, dis) can be 
uniquely determined in terms of the known constants 
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given by the left-hand side of Eas. ((36)) and (|37|. 
Similar conditions are obtained for the second func- 
tion {ul '{p),Ve (p))- However, the two functions 
{ug{p),Vg{p)) obtained in this way for all values of p are 
not properly normalized in the continuum and are not 
orthogonal to each other. In this case, the orthonormal- 
ization condition (fTOl) reads: 



(38) 



dpp u^;(pX(p) + v^,{p)v^{p) 



^ [cxicx'i + dxidx'i] S{k^l'^ - k^l^ ) 

[cx3Cx'3 + dxsdx'i] —Ty- 5{kf'' - kf^ ) 
k \ 



where A, A' = (1,2). 



APPENDIX B: REGULARIZATION 
PROCEDURE FOR THE SELF-CONSISTENT 
GAP EQUATION 

It is well known that the self-consistent condition ([2]) 
for the gap parameter A(r) diverges in the ultraviolet in 
the case of a contact inter-particle potential with coupling 
constant vq and has to be regularized accordingly. 

In the homogeneous case with a uniform gap parameter 
Aq, this rcgularization is readily achieved by expressing 
the bare coupling constant vq that enters Eq.Q in terms 
of the scattering length a^? of the two-body problem, via 
the relation 



1 

vq 



(27r)3 k2 



(39) 



Here, fco is an ultraviolet cutoff which is eventually let 
— > oo while t^o by keeping ap sX the desired value. 

This simple rcgularization, however, cannot be ex- 
ploited when the gap parameter A(r) has a spatial de- 
pendence occurring, for instance, in the presence of an 
isolated vortex as considered in the present paper, or, 
more generally, in the presence of a scalar trapping po- 
tential K,xt(r) or of an effective vector potential AjV), 



the latter arising when the trap is set into rotation 11 1 
or artificial gauge potentials are applied to neutral atoms 
[39} . In all these cases, a new strategy is required. 

A number of procedures have already been devised to 
implement a consistent rcgularization scheme for inhomo- 
geneous situations, ranging from the simple introduction 
of an energy cutoff, to relying on the pseudo-potential 
method to regularize the anomalous density in real space 
(40| . and to a combination of an energy cutoff with a 
local-density approximation [isj (which has then be sub- 
ject to improvements j20l|). 

In this Appendix, we introduce a procedure to regu- 
larize the gap equation ^ under generic inhomogeneous 
situations, which combines the introduction of an energy 



cutoff Ec as done in subsection II-A for the explicit nu- 
merical solution of the BdG equations for eigenvalues 
up to the value (Ec — p), with the derivation of the Gross- 
Pitaevskii equation for composite bosons that was done 
in Ref. |2l| in the BEG limit starting from the BdG equa- 
tions in terms of the small quantity A(r)/|/x|. Our rcg- 
ularization procedure for the gap equation ([2]), however, 
is not limited to the BEG limit but holds instead for any 
coupling throughout the BCS-BEC crossover, since in the 
present context it is the ratio A{r) /{Ec — p) to play the 
role of the small quantity that allows for the identification 
of the terms to be retained in the final expression. 

Schematically, our rcgularization procedure of the gap 
equation ^ is based on the following steps: 

(i) We consider a wave-vector cutoff fc^ such that Ec = 
kc/{2m) is the energy cutoff introduced in subsection II- 
A. We take Ec » Ep where Ep = kp/ (2m) is the Fermi 
energy associated with the mean density uq = A;|,/(37r^). 
While kc will be kept finite in the calculation, the ultra- 
violet cutoff fco ^ kc entering Eq. (p9|) will eventually be 
taken to diverge; 

(ii) We split the in Eq.(I2|) in two parts, with < 
Ec — p and > Ec ~ p, respectively. While in the 
first part the wave functions {u^(t),v^{t)) are explic- 
itly calculated numerically, the second part is treated 
within a local-density approximation as specified below. 
In addition, the Fermi function frif-v) can be dropped 
from this second part for all practical purposes, because 
Ec/ [ksT) 3> 1 even when ksT is of the order of Ep; 

(iii) Following Ref. jl8j. we rewrite the integral on the 
right-hand side of Eg. ([55]) as follows: 



ho 



dk m 



(27r)3 k2 



n{kc) 



dk 



1 



(27r)3 k2/m-2^ 



(40) 



that defines the quantity TZ{kc). A simple calculation 
then yields: 



2^2 



nikc) 



/ kc — %/27n/i \ 
\ kc-j'^/2m|J, J 



(m>o) 



kc + 2m\i 



T7 — arctan , , 



(m < 0) . 
(41) 

Through the above steps, the self-consistent condition 
^ for the gap parameter becomes: 



n{kc) 



1 



dk 

(27r)3 k2/m-2^ 



A(r) 



V 

5Z ^^'^(i")w..(r)* . 



(42) 



Here, like in Ref.[18[, the last term within brackets on the 
left-hand side of Eq. ip^ can be used to regularize the 
with e,y > Ec — p on the right-hand side. What is novel 
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of the present approach, however, is the way this high- 
energy sum is dealt with, by drawing connections with 
the derivation of the Gross-Pitaevskii equation for com- 
posite bosons in the BEG hmit that was done in Ref. [2lj 
starting from the fuU BdG equations. 

To this end, we refer directly to Eq.(13) of Ref. pll and 
write (by also keeping the same notation of Ref. |21| ) : 



e„>Bc- 

E 



y"drig(r,ri)* A(ri) 



(43) 



/ 



dridradra i?(r, ri, ra, rg)* A(ri) A(r2)* A(r3) 



where Q(r,ri) and i?(r, ri, r2, ra) are defined by 
Eqs.(14) and (15) of Ref. 2l|, in the order, but are here 
considered with the provision that all h-integrals that en- 
ter those expressions through the Fourier representation 
of the non-interacting Green's function Qq therein are re- 
stricted by |k| > kc- In addition, the local-density con- 
ditionu -> ^(r) = /i — Vc^t{i^) is adopted here like in 
Ref.|2l|. to take into account the possible presence of a 
trapping potential. 

Following further Eq.(16) of Ref. [2l| . we approximate 
for a sufficiently slowly varying gap parameter A(r): 



y"driQ(r,ri)*A(ri) = 



ao(r)' 



where now 
ao(r) 



A(r) 

(44) 



|k|>,,(2^)3 kVm- 
dk 1 



2^-h2V;xt(r) 



|k|>,, (2^)3 kVm - 2^^ 
2l^oxt(r) ' 



|k|>fcj2^)' (k2/m-2^)' 



(45) 



and 



|k|>fc, 
1 



Am 



2m 



6m 



(46) 



Note that, once these expressions arc used in Eq. p2)) . the 
first term on the right-hand side of Eg. ipS)) cancels the 
last term within brackets on the left-hand side of Eq. ([^^ . 
By a similar token we obtain: 



j dvxdY2dvs R{r, ri , r2, ra) 



dk 



Introducing at this point the notation 

dk \2m J 

|k|>fc,(2^ 



< 2m 



(47) 



(48) 



the expression p3)) can be written compactly as follows: 



dk 



k|>fej2^)3 k2/m-2/i 



-214xt(r)ilo2(fcc) 



■Io3(fcc)|A(r)pA(r) 



V^A(r) 

Am 



A(r) 



(49) 



Here, the first term on the right-hand side which is linear 
in A(r) was already introduced in Ref. 18|, while the ad- 
dition of the second term on the right-hand side which is 
cubic in A(r) was already considered in Ref. (20t . What 
comes naturally from the present derivation is the fur- 
ther introduction of the third term on the right-hand 
side which emphasizes the spatial variations of A(r). 

Entering the approximate expression (j49p into the 
right-hand side of Eq. (1421) yields eventually the regular- 
ized gap equation we were looking for: 



Airaf 



^Io2{kc) - ^Iisikc 



Am 



2Kxt(r)ilo2(fcc) + ilo3(fcc)|A(r)|2 ^A(r) 
u,{r)vAr)* [I - 2fF{e,)] . 



(50) 



Note that the expression on the right-hand side, which 
results from an explicit numerical integration of the BdG 
equations, acts as a source term on the non-linear differ- 
ential equation for A(r) given by the left-hand side. 

But for the source term on its right-hand side, Ea. ([5n)) 
resembles the Gross-Pitaevskii equation with suitable co- 
efficients, and actually reduces to it in the BEG limit 
when ^ (< 0) is the largest energy scale in the problem, 
such that Ec (and thus kc) can be taken to vanish for 
all practical purposes. In particular, in this limit one 
obtains for the integrals entering Eq. (|50| the values: 



Hike) 

2^13 (^c' 



m^ap 



Anap Stt 



fiB , 2:02 (fcc) 



m^ap 
2tt 



Stt 



3 3 
map 

An 



(51) 



where = 2/i + (ma^)^^ is the chemical potential for 
composite bosons. With the rescaling given in Ref. (2l| . 
between the gap function A(r) and the condensate wave 
function for composite bosons, one recovers in this way 
from Eq. ([50| the Gross-Pitaevskii equation given by 
Eq.(20) of Ref.[2l|. 

In practice, the inclusion of successively more terms on 
the right-hand side of Eq. p9|) [from the linear (A) term, 
to the linear plus cubic (A + A^) terms, and finally to 
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LJJ 




case, the smaller value = 3Ep of the cutoff with re- 
spect to the value Ec = 9Ep (or Ec = ISEp, depending 
on the coupling) needed to the full calculation for achiev- 
ing a stable configuration, yields a reduction of the total 
computational time by a factor of five or more. 

Specifically, one sees from Fig[T4] that inclusion of all 
terms on the right-hand side of Eq. p9)) [namely, the lin- 
ear plus cubic plus Laplacian (A -f- A'^ -|- V^) terms] leads 
for all couplings to quite a good agreement with the full 
calculation, and not only in the asymptotic (bulk) re- 
gion but also near the center of the vortex where A(p) is 
strongly depressed. In contrast, the approximation that 
includes only the first two terms on the right-hand side 
of Eq. p^ [namely, the linear plus cubic (A-|- A'^) terms] 
reproduces the bulk value Aq but leads to (even sizable) 
deviations from the full calculation near the center of the 
vortex. Finally, the approximation that includes only the 
first term on the right-hand side of Eq. ([1^1) [namely, the 
linear (A) term] progressively deviates for all values of 
p from the full calculation when approaching the BEC 
limit, where it is not able to recover the correct bulk 



value Aq with the required accuracy [41 



APPENDIX C: ASYMPTOTIC FORM OF THE 
NUMBER AND CURRENT DENSITIES 



FIG. 14. A(p) (in units of Ef) vs p (in units of kp^) ob- 
tained by various approximations at zero temperature for: 



(a) {kpap) 



-1; (b) {kpap) 



0; (c) (kpap)' 



Comparison is made between the "best" calculation with a 
large value of Ec [that equals 9Ep in panels (a) and (b), and 
18Ep in panel (c)], and less sophisticated calculations all with 
the smaller value Ec = SEp which include, respectively, the 
linear, the linear plus cubic, and the linear plus cubic plus 
Laplacian terms on the right-hand side of Eq. 1)49^ . 



the linear plus cubic plus Laplacian (A + A'^ + V^) terms] 
enables one to decrease the total computational time at 
any coupling by decreasing the value of the cutoff Ec up 
to which the eigenfunctions of the BdG equations have 
to be explicitly calculated. This can be achieved without 
loosing accuracy in the shape of A(r) as well as of other 
physical quantities (see also Appendix C). 

As an example, we consider again the problem of an 
isolated vortex in an otherwise infinite superfluid, which 
is the main concern of the present paper. For this case, 
the profile of A(p) vs p at zero temperature for three dif- 
ferent couplings across the BCS-BEC crossover is shown 
in FigfT4l where alternative numerical approximations 
(including the linear, linear plus cubic, and linear plus 
cubic plus Laplacian terms), which all adopt a common 
and rather small value of the cutoff Ec, are compared 
with the full calculation of the BdG equations where the 
value of the cutoff Ec is taken considerably larger. In this 



Besides the gap A(r), other relevant physical quanti- 
ties obtained by solving the BdG equations ^ are the 
number n(r) and current j(r) densities. They are given, 
respectively, by the expressions: 

n(r) = 2^ [fp{e^)\u^{r)\' + (1 - fp{e^)) \v4r)\'] (52) 

j(r) = — V{/F(e.) [u,(r)*Vn.(r) - (Vu4r)*)u.(r)] 
im ^-^ 

+ (1 - .^(e,.)) [v,{r)^v,{v)* - (Vw4r))i;,(r)*]} . (53) 

Only positive eigenvalues can be considered for the sums 
in Eqs.dll) and (p]). 

In order to calculate the expressions ([52]) and ([53| in 
an efficient way, and consistently with what was done 
in Appendix B for the gap equation also the in 
Eas. ((52|) and ([53]) is split into two parts, with < Ec — p. 
and £i, > Ec — p. While in the first part with < Ec — p 
one uses explicitly the wave functions (u^(r), 7j,y(r)) ob- 
tained by solving the BdG equations, for both quan- 
tities 7i(r) and j(r) the second (asymptotic) part with 
£i, > Ec — p is treated again within a local-density ap- 
proximation as follows. 

Similarly to what was done in Appendix B, we adapt 



to the present situation the treatment made in Ref.[21 



where expressions for n(r) and j(r) consistent with the 
Gross-Pitaevskii equation were recovered in the BEC 
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FIG. 15. (a) Number density n(p) (normalized to its bulk 
value no) and (b) current density j(p) (normalized to its max- 
imum value at Rv) vs p (in units oi kp^) obtained at zero tem- 
perature and unitarity for an isolated vortex. Calculations 
corresponding to two values of the cutoff Ec are compared 
with each other. 
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FIG. 16. The values of the radius i?v of the vortex (dots) 
are shown versus the temperature T (in units of the respec- 
tive critical temperature Tc) for four different couplings across 
unitarity. These values are then fitted by the mean-field-like 
expression Rv{T) oc {Tc — T)~^^'^ over the whole temperature 
range down to T = (dashed lines). 



limit. Wc thus obtain for the asymptotic parts: 

nasy„.(r)= b.(r)|'= 2lo2(M|A(r)r (54) 

and 

jasym(r) = — V K (r) V (r) * - (Vw^ (r))w^ (r) *] 

V 

X [A(r)*VA(r) - A(r)VA(r)*] (55) 

with the definition pS)) for the integrals over the wave 
vector k with |k| > k^. 

The above expressions hold for any coupling. In partic- 
ular, in the BEC limit, whereby the integrals Io2(fcc) and 
2^03 (fee) reduce to the values ((5T|) . the expressions ((54)) and 
(fF5|) reduce to Eqs.(21) and (22) of Ref.jlH, in t he order, 
once the proper rescaling $(r) = y'm'^aF / {^t^) A(r) be- 
tween the gap function and the condensate wave function 
$(r) for composite bosons is performed. 



The above expressions can again be applied to the 
problem of an isolated vortex in an otherwise infinite su- 
perfluid. In particular, in Fig llSI we show the results of 
our calculation for the number density n[p) and the cur- 
rent density j(p) at a distance p from the center of the 
vortex, when different values of the cutoff are used. 
The calculation is done at zero temperature and unitar- 
ity. Once again we verify that, with the complete regu- 
larization procedure for the gap equation introduced in 
Appendix B and here extended to the density and cur- 
rent, rather small values of Ec are sufficient in practice to 
obtain results in quite good agreement with our "best" 
calculation where a large value of Ec is used. In addi- 
tion, the results of Fig llSI may serve also implicitly to 
verify that the approximate expressions (|54p and (|55p . 
respectively for the asymptotic density and current, were 
obtained in a physically sound manner. 

Note finally that, by this kind of plots, the radius i?v 
of the vortex can be identified as corresponding to the 
maximum value of the current. In turn, the value of 
i?v fixes a length scale which is relevant to reckon the 
value of i?out J that was introduced in subsection II- A and 
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used in Appendix A to enforce the boundary conditions 
on the radial wave functions (typical values of the ratio 
-Rout / taken in the calculations range from about 50 at 
low temperature to about 10 close to for the couplings 
we have explored). 

The values of obtained in this way are reported 
in Fig |16l for four couplings across unitarity versus the 
temperature T (in units of the respective critical tem- 
perature Tc). These values are then fitted by the mean- 
field-like expression kpR^ = B {1 — T/Tc)'^^"^ , obtaining 
the values B = (17.72,4.26, 1.41, 1.08) for the couplings 
{kpap)"^ (-2.0, -1.0, 0.0, -1-1.0), in the order. On the 
average, these values for the prc-factor B are larger by 
•\/2 than the values for the prc-factor A entering the cor- 
responding expression kp£^{T) ^ A{1 — T/Tc)^^^^ that 
were reported in Section IV, thus confirming our conclu- 
sion made also in Section IV that a single length scale 
can be extracted from the BdG equations. 
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